Practical - Week 2
Florencia Grattarola
(Department of Spatial Sciences)
2022-10-03
Gridding biodiversity data, exploring patterns at different resolutions, and making pretty maps.
Examples
For all geocomputationgs today we will use the sf package.
For all the visualisations today we will use the ggplot package. It’s part of the tidyverse.
Photo of Capreolus capreolus observed in Czech Republic by romanvrbicek licensed CC-BY-NC, via iNaturalist
Steps:
POINTS)Get the occurrence records
# A tibble: 3,050 × 174
key scien…¹ decim…² decim…³ issues datas…⁴ publi…⁵ insta…⁶ publi…⁷ proto…⁸
<chr> <chr> <dbl> <dbl> <chr> <chr> <chr> <chr> <chr> <chr>
1 34557… Myodes… 49.8 16.0 cdrou… 50c950… 28eb1a… 997448… US DWC_AR…
2 34558… Sus sc… 49.2 16.5 cdrou… 50c950… 28eb1a… 997448… US DWC_AR…
3 34558… Myotis… 50.2 16.8 cdrou… 50c950… 28eb1a… 997448… US DWC_AR…
4 34558… Lepus … 50.3 15.3 cdrou… 50c950… 28eb1a… 997448… US DWC_AR…
5 34558… Capreo… 50.3 15.3 cdrou… 50c950… 28eb1a… 997448… US DWC_AR…
6 34561… Cervus… 50.3 14.6 cdrou… 50c950… 28eb1a… 997448… US DWC_AR…
7 34562… Capreo… 49.2 17.4 cdrou… 50c950… 28eb1a… 997448… US DWC_AR…
8 34564… Dama d… 49.2 16.5 cdrou… 50c950… 28eb1a… 997448… US DWC_AR…
9 34564… Lepus … 50.3 15.3 cdrou… 50c950… 28eb1a… 997448… US DWC_AR…
10 34565… Myocas… 50.0 14.7 cdrou… 50c950… 28eb1a… 997448… US DWC_AR…
# … with 3,040 more rows, 164 more variables: lastCrawled <chr>,
# lastParsed <chr>, crawlId <int>, hostingOrganizationKey <chr>,
# basisOfRecord <chr>, occurrenceStatus <chr>, taxonKey <int>,
# kingdomKey <int>, phylumKey <int>, classKey <int>, orderKey <int>,
# familyKey <int>, genusKey <int>, speciesKey <int>, acceptedTaxonKey <int>,
# acceptedScientificName <chr>, kingdom <chr>, phylum <chr>, order <chr>,
# family <chr>, genus <chr>, species <chr>, genericName <chr>, …
Let’s transform our data table into an sf object using st_as_sf()
Let’s transform our data table into an sf object using st_as_sf()
mammalsCZ %>%
filter(!is.na(decimalLongitude) & !is.na(decimalLatitude)) %>% # filter records without coordinates
st_as_sf(coords=c('decimalLongitude', 'decimalLatitude'),
crs=4326)Simple feature collection with 3049 features and 6 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 12.20682 ymin: 48.60141 xmax: 18.79605 ymax: 50.99892
Geodetic CRS: WGS 84
# A tibble: 3,049 × 7
species order event…¹ state…² count…³ datas…⁴ geometry
* <chr> <chr> <chr> <chr> <chr> <chr> <POINT [°]>
1 Myodes glare… Rode… 2022-0… Pardub… CZ iNatur… (15.95647 49.80115)
2 Sus scrofa Arti… 2022-0… Jihomo… CZ iNatur… (16.51939 49.19975)
3 Myotis myotis Chir… 2022-0… Pardub… CZ iNatur… (16.84201 50.18355)
4 Lepus europa… Lago… 2022-0… Středo… CZ iNatur… (15.33439 50.28446)
5 Capreolus ca… Arti… 2022-0… Středo… CZ iNatur… (15.33636 50.28574)
6 Cervus nippon Arti… 2022-0… Středo… CZ iNatur… (14.58019 50.3448)
7 Capreolus ca… Arti… 2022-0… Zlínský CZ iNatur… (17.36916 49.23076)
8 Dama dama Arti… 2022-0… Jihomo… CZ iNatur… (16.53508 49.2115)
9 Lepus europa… Lago… 2022-0… Králov… CZ iNatur… (15.2629 50.30738)
10 Myocastor co… Rode… 2022-0… Středo… CZ iNatur… (14.653 49.99369)
# … with 3,039 more rows, and abbreviated variable names ¹eventDate,
# ²stateProvince, ³countryCode, ⁴datasetName
Note that crs=4326 is related to the EPSG:4326 CRS: WGS 84
Let’s transform our data table into an sf object using st_as_sf()
Get the administrative borders
Reading layer `CZE_adm0' from data source
`/Users/flograttarola/Documents/GitHub/Spatial-Ecology-and-Macroecology/Practical_classes/Week2_gridding_and plotting/data/CZE_adm0.shp'
using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 70 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 12.08586 ymin: 48.54084 xmax: 18.86253 ymax: 51.05438
Geodetic CRS: WGS 84
Simple feature collection with 1 feature and 70 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 12.08586 ymin: 48.54084 xmax: 18.86253 ymax: 51.05438
Geodetic CRS: WGS 84
GADMID ISO NAME_ENGLI NAME_ISO NAME_FAO NAME_LOCAL
1 59 CZE Czech Republic CZECH REPUBLIC Czech Republic Ceská republika
NAME_OBSOL NAME_VARIA NAME_NONLA NAME_FRENC NAME_SPANI
1 Moravia|Bohemia Cesko <NA> République Tchèque República Checa
NAME_RUSSI NAME_ARABI NAME_CHINE WASPARTOF CONTAINS
1 ????? ????????? ???????? ????? Czechoslovakia <NA>
SOVEREIGN ISO2 WWW FIPS ISON VALIDFR VALIDTO AndyID POP2000 SQKM
1 Czech Republic CZ <NA> EZ 203 19930101 Present 65 10271830 78495.16
POPSQKM UNREGION1 UNREGION2 DEVELOPING CIS Transition OECD
1 130.8594 Eastern Europe Europe 2 0 0 1
WBREGION WBINCOME WBDEBT WBOTHER CEEAC CEMAC
1 Europe & Central Asia Upper middle income Less indebted <NA> 0 0
CEPLG COMESA EAC ECOWAS IGAD IOC MRU SACU UEMOA UMA PALOP PARTA CACM EurAsEC
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
Agadir SAARC ASEAN NAFTA GCC CSN CARICOM EU CAN ACP Landlocked AOSIS SIDS
1 0 0 0 0 0 0 0 1 0 0 1 0 0
Islands LDC Shape_Leng Shape_Area geometry
1 0 0 26.21557 9.813959 MULTIPOLYGON (((13.17648 49...
Get the standard grids from the Czech Republic
Reading layer `KvadratyCR_JTSK' from data source
`/Users/flograttarola/Documents/GitHub/Spatial-Ecology-and-Macroecology/Practical_classes/Week2_gridding_and plotting/data/KvadratyCR_JTSK.shp'
using driver `ESRI Shapefile'
Simple feature collection with 678 features and 10 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: -912065.7 ymin: -1229672 xmax: -421278.2 ymax: -928997.8
Projected CRS: S-JTSK / Krovak East North
Simple feature collection with 678 features and 10 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: -912065.7 ymin: -1229672 xmax: -421278.2 ymax: -928997.8
Projected CRS: S-JTSK / Krovak East North
First 10 features:
OBJECTID ENTITY AREA PERIMETE KVADRAT X Y Shape_Leng
1 1 3 130.032 45.626 4951 -740241.8 -935317.5 45626.41
2 2 5 130.036 45.627 4952 -728665.0 -936923.5 45627.17
3 3 6 130.303 45.675 5051 -741782.9 -946335.6 45675.21
4 4 8 130.040 45.628 4953 -717084.4 -938503.6 45627.96
5 5 9 130.307 45.676 5052 -730181.5 -947945.1 45675.88
6 6 12 130.311 45.677 5053 -718576.3 -949528.8 45676.59
7 7 14 130.050 45.630 4955 -693912.1 -941586.4 45629.61
8 8 17 130.055 45.630 4956 -682320.5 -943089.1 45630.48
9 9 18 130.319 45.678 5055 -695354.9 -952618.5 45678.09
10 10 19 130.060 45.631 4957 -670725.5 -944566.0 45631.36
Shape_Area ID geometry
1 130031546 4951 POLYGON ((-745252.4 -928997...
2 130035855 4952 POLYGON ((-733689.7 -930614...
3 130302745 5051 POLYGON ((-746805.7 -940014...
4 130040349 4953 POLYGON ((-722123.3 -932206...
5 130306551 5052 POLYGON ((-735218.5 -941635...
6 130310574 5053 POLYGON ((-723627.4 -943229...
7 130049749 4955 POLYGON ((-698979.2 -935311...
8 130054715 4956 POLYGON ((-687401.8 -936825...
9 130319128 5055 POLYGON ((-700434.2 -946342...
10 130059749 4957 POLYGON ((-675820.7 -938313...
Check the layer’s Coordinate Reference System (CRS)
Coordinate Reference System:
User input: WGS 84
wkt:
GEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["latitude",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["longitude",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]]
Coordinate Reference System:
User input: S-JTSK / Krovak East North
wkt:
PROJCRS["S-JTSK / Krovak East North",
BASEGEOGCRS["S-JTSK",
DATUM["System of the Unified Trigonometrical Cadastral Network",
ELLIPSOID["Bessel 1841",6377397.155,299.1528128,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4156]],
CONVERSION["Krovak East North (Greenwich)",
METHOD["Krovak (North Orientated)",
ID["EPSG",1041]],
PARAMETER["Latitude of projection centre",49.5,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8811]],
PARAMETER["Longitude of origin",24.8333333333333,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8833]],
PARAMETER["Co-latitude of cone axis",30.2881397527778,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",1036]],
PARAMETER["Latitude of pseudo standard parallel",78.5,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8818]],
PARAMETER["Scale factor on pseudo standard parallel",0.9999,
SCALEUNIT["unity",1],
ID["EPSG",8819]],
PARAMETER["False easting",0,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["easting (X)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["northing (Y)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["GIS."],
AREA["Czechia; Slovakia."],
BBOX[47.73,12.09,51.06,22.56]],
ID["EPSG",5514]]
We will transform mammalsCZ CRS to S-JTSK / Krovak East North using st_transform()
Simple feature collection with 3049 features and 6 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -896937 ymin: -1224642 xmax: -435147.4 ymax: -943016.8
Projected CRS: S-JTSK / Krovak East North
# A tibble: 3,049 × 7
species order event…¹ state…² count…³ datas…⁴ geometry
* <chr> <chr> <chr> <chr> <chr> <chr> <POINT [m]>
1 Myodes glare… Rode… 2022-0… Pardub… CZ iNatur… (-637524.1 -1088504)
2 Sus scrofa Arti… 2022-0… Jihomo… CZ iNatur… (-604525.8 -1159543)
3 Myotis myotis Chir… 2022-0… Pardub… CZ iNatur… (-569727 -1053229)
4 Lepus europa… Lago… 2022-0… Středo… CZ iNatur… (-675301.1 -1029784)
5 Capreolus ca… Arti… 2022-0… Středo… CZ iNatur… (-675143.9 -1029662)
6 Cervus nippon Arti… 2022-0… Středo… CZ iNatur… (-727699 -1016188)
7 Capreolus ca… Arti… 2022-0… Zlínský CZ iNatur… (-542595 -1162493)
8 Dama dama Arti… 2022-0… Jihomo… CZ iNatur… (-603247.5 -1158367)
9 Lepus europa… Lago… 2022-0… Králov… CZ iNatur… (-680037.4 -1026620)
10 Myocastor co… Rode… 2022-0… Středo… CZ iNatur… (-727767.4 -1055585)
# … with 3,039 more rows, and abbreviated variable names ¹eventDate,
# ²stateProvince, ³countryCode, ⁴datasetName
Let’s plot the sf objects
To calculate number of records (N) and number of species per grid-cell (SR), we will use st_join()
To calculate number of records (N) and number of species per grid-cell (SR), we will use st_join().
Simple feature collection with 3350 features and 16 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: -912065.7 ymin: -1229672 xmax: -421278.2 ymax: -928997.8
Projected CRS: S-JTSK / Krovak East North
First 10 features:
OBJECTID ENTITY AREA PERIMETE KVADRAT X Y Shape_Leng
1 1 3 130.032 45.626 4951 -740241.8 -935317.5 45626.41
2 2 5 130.036 45.627 4952 -728665.0 -936923.5 45627.17
3 3 6 130.303 45.675 5051 -741782.9 -946335.6 45675.21
4 4 8 130.040 45.628 4953 -717084.4 -938503.6 45627.96
5 5 9 130.307 45.676 5052 -730181.5 -947945.1 45675.88
5.1 5 9 130.307 45.676 5052 -730181.5 -947945.1 45675.88
5.2 5 9 130.307 45.676 5052 -730181.5 -947945.1 45675.88
5.3 5 9 130.307 45.676 5052 -730181.5 -947945.1 45675.88
5.4 5 9 130.307 45.676 5052 -730181.5 -947945.1 45675.88
5.5 5 9 130.307 45.676 5052 -730181.5 -947945.1 45675.88
Shape_Area ID species order eventDate
1 130031546 4951 <NA> <NA> <NA>
2 130035855 4952 <NA> <NA> <NA>
3 130302745 5051 <NA> <NA> <NA>
4 130040349 4953 <NA> <NA> <NA>
5 130306551 5052 Capreolus capreolus Artiodactyla 2020-06-23T00:00:00
5.1 130306551 5052 Capreolus capreolus Artiodactyla 2020-06-29T00:00:00
5.2 130306551 5052 Capreolus capreolus Artiodactyla 2020-06-27T00:00:00
5.3 130306551 5052 Lepus europaeus Lagomorpha 2020-06-29T00:00:00
5.4 130306551 5052 Capreolus capreolus Artiodactyla 2020-06-24T00:00:00
5.5 130306551 5052 Lepus europaeus Lagomorpha 2020-06-26T00:00:00
stateProvince countryCode datasetName geometry
1 <NA> <NA> <NA> POLYGON ((-745252.4 -928997...
2 <NA> <NA> <NA> POLYGON ((-733689.7 -930614...
3 <NA> <NA> <NA> POLYGON ((-746805.7 -940014...
4 <NA> <NA> <NA> POLYGON ((-722123.3 -932206...
5 <NA> CZ <NA> POLYGON ((-735218.5 -941635...
5.1 <NA> CZ <NA> POLYGON ((-735218.5 -941635...
5.2 <NA> CZ <NA> POLYGON ((-735218.5 -941635...
5.3 <NA> CZ <NA> POLYGON ((-735218.5 -941635...
5.4 <NA> CZ <NA> POLYGON ((-735218.5 -941635...
5.5 <NA> CZ <NA> POLYGON ((-735218.5 -941635...
Now, we need to summarise the data per grid-cell.
Luckily, Tidyverse methods also work for sf objects :)
We will do it with group_by() and summarise().
First, we group by grid-cells and then we count values per grid.
Let’s calculate number of records (N) and number of species per grid-cell (SR)
Simple feature collection with 678 features and 3 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: -912065.7 ymin: -1229672 xmax: -421278.2 ymax: -928997.8
Projected CRS: S-JTSK / Krovak East North
# A tibble: 678 × 4
OBJECTID N SR geometry
<int> <int> <int> <POLYGON [m]>
1 1 1 1 ((-745252.4 -928997.8, -733689.7 -930614.9, -733689.7 -…
2 2 1 1 ((-733689.7 -930614.9, -722123.3 -932206.2, -722123.3 -…
3 3 1 1 ((-746805.7 -940014.3, -735218.5 -941635, -735218.5 -94…
4 4 1 1 ((-722123.3 -932206.2, -710553.1 -933771.6, -710553.1 -…
5 5 27 3 ((-735218.5 -941635, -723627.4 -943229.9, -723627.4 -94…
6 6 4 3 ((-723627.4 -943229.9, -712032.7 -944798.9, -712032.7 -…
7 7 1 1 ((-698979.2 -935311.3, -687401.8 -936825.2, -687401.8 -…
8 8 1 1 ((-687401.8 -936825.2, -675820.7 -938313.3, -675820.7 -…
9 9 1 1 ((-700434.2 -946342, -688832.2 -947859.3, -688832.2 -94…
10 10 1 1 ((-675820.7 -938313.3, -664236.1 -939775.7, -664236.1 -…
# … with 668 more rows
Let’s store the output into a new object CZ_mammals_grids.
Finally, let’s plot this.
We will do it using geom_sf() a ggplot function to visualise sf objects.
Where are the nice colours?
We need to indicate which column from the object should the grids be filled with.
Let’s get a nicer color scale.
Bare in mind that we see colors differently. Thus, it’s important to consider colorblind safe palettes
Check out The R Graph Gallery for many cool graphic options with ggplot
ggplot() +
geom_sf(data=mammalsCZ_grids, aes(fill=SR)) +
scale_fill_fermenter(palette ='YlOrBr', n.breaks=9, direction = 1)+
geom_sf(data=CZ_borders, fill=NA) +
coord_sf(crs=4326) +
theme_bw()Now, here’s a better plot :)
The hotspots are located in cities, how can they be the richest areas for mammals?
Let’s have a look at the sampling effort (N: number of records) and compare both layers.
What can you say about the hotspots of species richness we found?
POINTS)
We are done! No it’s your turn :)
Photo of Emys orbicularis observed in France by Martin Costechareire licensed CC-BY-NC, via iNaturalist
Steps:
POLYGONS)Read testudines’ IUCN range maps (POLYGONS)
Reading layer `data_0' from data source
`/Users/flograttarola/Documents/GitHub/Spatial-Ecology-and-Macroecology/Practical_classes/Week2_gridding_and plotting/data/testudines/data_0.shp'
using driver `ESRI Shapefile'
Simple feature collection with 169 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -179.999 ymin: -47.55406 xmax: 179.999 ymax: 61.50515
Geodetic CRS: WGS 84
Polygons were downloaded from IUCN Spatial Data Download’s page.
What do the data look like?
Simple feature collection with 169 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -179.999 ymin: -47.55406 xmax: 179.999 ymax: 61.50515
Geodetic CRS: WGS 84
First 10 features:
ASSESSMENT ID_NO BINOMIAL PRESENCE ORIGIN SEASONAL
1 495907 10041 Heosemys annandalii 1 1 1
2 499158 10825 Indotestudo forstenii 1 1 1
3 499618 10950 Pangshura sylhetensis 1 1 1
4 505402 12124 Lissemys scutata 1 1 1
5 507698 12695 Malaclemys terrapin 1 1 1
6 507698 12695 Malaclemys terrapin 1 5 1
7 508210 12696 Malacochersus tornieri 1 1 1
8 508518 12775 Manouria impressa 1 1 1
9 511526 13038 Melanochelys tricarinata 1 1 1
10 511526 13038 Melanochelys tricarinata 6 1 1
COMPILER YEAR CITATION
1 Anders Rhodin 2020 Chelonian Research Foundation
2 Anders Rhodin 2020 Chelonian Research Foundation
3 Anders Rhodin 2020 Chelonian Research Foundation
4 Anders Rhodin 2020 Chelonian Research Foundation
5 Anders Rhodin 2018 Chelonian Research Foundation
6 Anders Rhodin 2018 Chelonian Research Foundation
7 Anders Rhodin 2018 Chelonian Research Foundation
8 Anders Rhodin 2020 Chelonian Research Foundation
9 Anders Rhodin 2018 Chelonian Research Foundation
10 Anders Rhodin 2018 Chelonian Research Foundation
LEGEND SUBSPECIES SUBPOP DIST_COMM ISLAND
1 Extant (resident) <NA> <NA> <NA> <NA>
2 Extant (resident) <NA> <NA> <NA> <NA>
3 Extant (resident) <NA> <NA> <NA> <NA>
4 Extant (resident) <NA> <NA> <NA> <NA>
5 Extant (resident) <NA> <NA> <NA> <NA>
6 Extant & Origin Uncertain (resident) <NA> <NA> <NA> <NA>
7 Extant (resident) <NA> <NA> <NA> <NA>
8 Extant (resident) <NA> <NA> <NA> <NA>
9 Extant (resident) <NA> <NA> <NA> <NA>
10 Presence Uncertain <NA> <NA> <NA> <NA>
TAX_COMM geometry
1 <NA> MULTIPOLYGON (((102.6212 5....
2 <NA> MULTIPOLYGON (((120.8806 1....
3 <NA> MULTIPOLYGON (((92.29536 22...
4 <NA> MULTIPOLYGON (((95.29987 15...
5 <NA> MULTIPOLYGON (((-81.45417 2...
6 <NA> MULTIPOLYGON (((-64.7296 32...
7 <NA> MULTIPOLYGON (((32.93152 -9...
8 <NA> MULTIPOLYGON (((101.684 5.0...
9 <NA> MULTIPOLYGON (((85.58097 25...
10 <NA> MULTIPOLYGON (((92 21.59583...
Now, let’s get a world map. For this we will use the package rnaturalearth.
We will download a shapefile of the world at medium scale and we will combine all countries into one unique polygon.
Let’s see how it looks
Before working with this map, we need to project the layer: from lat/lon to equal area projection.
Earth is not flat :) Projections help us represent the two-dimensional curved surface of the globe into 2D space. There are many ways to do this (map projections).
Equal-area maps preserve area measure, generally distorting shapes in order to do that.
We will choose Equal Earth (EPSG:8857).
Let’s transform the data (both the world and the testudines’ layers)
Now, let’s create 1-degree grid-cells for the entire world.
We will do this using the function st_make_grid()
For the cellsize we will chose 1 degree, which is ~100km (=100,000m)
We have the grid, we are ready to calculate richness.
To make things faster, let’s get a smaller dataset of Testudines
This can take a while
Let’s plot 10 species as an example
To calculate number of records (N) and number of species per grid-cell (SR), we will use st_join().
This will take a while
And now, we plot!
Here’s the results :)
Simple feature collection with 20146 features and 3 fields
Geometry type: GEOMETRY
Dimension: XY
Bounding box: xmin: -16920570 ymin: -8392928 xmax: 16920570 ymax: 8315130
Projected CRS: WGS 84 / Equal Earth Greenwich
# A tibble: 20,146 × 4
gridID N SR .
<int> <dbl> <int> <GEOMETRY [m]>
1 1 0 0 POINT (-16920565 -2061608)
2 2 2 2 MULTIPOLYGON (((-16903934 -2109412, -16904928 -2110323, -…
3 3 2 2 MULTIPOLYGON (((-16905520 -2108496, -16907015 -2103886, -…
4 4 2 2 POLYGON ((-16782581 -2193350, -16785819 -2191050, -167869…
5 5 3 3 POLYGON ((-16805308 -1830105, -16810425 -1830932, -168160…
6 6 2 2 POLYGON ((-16783058 -2419039, -16784544 -2418515, -167857…
7 7 2 2 MULTIPOLYGON (((-16779389 -2191507, -16775928 -2194312, -…
8 8 2 2 MULTIPOLYGON (((-16749600 -2289778, -16749039 -2290664, -…
9 9 2 2 POLYGON ((-16692191 -601861.4, -16691700 -601711.1, -1669…
10 10 2 2 POLYGON ((-16658328 -2432319, -16657130 -2432783, -166574…
# … with 20,136 more rows